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Abstract 



Presented in this paper is a new, efficient discrete-time domain system identification 
and model reduction method for large-scaled linear dynamic systems with multiple 
inputs. The method is based on a modification of the classical Eigensystem Real- 
ization Algorithm (ERA) and a simultaneous injection of multiple inputs, so called 
the Single-Composite-Input (SCI). Since the system response is sampled almost ex- 
clusively for the single representative input instead of the multiple individual inputs, 
this technique can significantly reduce the model construction time as well as the 
amount of the sampled data. For derivation of the new algorithm, a singular value 
decomposition is performed using a single set of output measurements that are not 
necessarily attributed to pulse inputs. Application to general computational fluid dy- 
namic systems and formulation of reduced-order aeroelastic models are also presented. 
Aerodynamic and aeroelastic systems modeled by the vortex lattice and CFL3D code 
are presented for a demonstration of the method. Based on the reduced computation 
time and excellent results obtained from the reduced-order models, the proposed al- 
gorithm shows a great potential as a linear system identification and model reduction 
tool for large-scaled systems subjected to multiple inputs. 



Nomenclature 



ABCD 



System matrices 
Aeroelastic system matrices 



A<n A d2 
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A s B s C s Structural system matrices 

b Reference length 

Cij Cross-correlation coefficient 

E Expected value 

K Covariance matrix denned in (48) 

L Dimension of original system 

M Number of time or frequency samples 

or Mach number 

mck Mass, damping, stiffness matrices 

p {Ri x 1) generalized coordinate vector 

q Dynamic pressure (= \pV 2 ) 

R Number of chosen singular modes or 

the dimension of realized model 

Ri Number of chosen KL modes 

s Laplace variable 

t Real time 

u (Ni x 1) input or generalized structural 

coordinate vector 

Ui i-th input or structural coordinate 

V Air speed 

x (L x 1) state or aerodynamic state vector 

X Fourier transform of x 

y (N Q x 1) output vector or (iVj x 1) 

generalzied aerodynamic force vector 

y" Pulse response due to i - th input 

<j>i KL mode 

$ KL modal matrix 

p Air density 
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t Reduced time (= ™) 

00 Frequency (rad/sec) 

uj c Maximum cut-off frequency 

Subscripts 

1 Input 

o Output 

R Reduced 

ref Reference 

s Structure 



1 Introduction 

Most modern dynamic models are constructed based on finite spatial discretizations 
of continuous systems, often resulting in a considerable number of degrees of freedom 
in the model. Consequently, for fast and efficient estimation of dynamic behavior 
as well as optimizations and closed-loop control designs, a model reduction must be 
accompanied. Essential requirements for reduced-order dynamic model are that the 
size of the system must not be too large, the model must be robust and have a good 
fidelity, it must be in the state-space, time domain formulation for implementation of 
active control systems and nonlinear time analysis, and finally, the reduction process 
itself must not be too expensive. 

There have been many model reduction methods available, but most of them 
require modifying the main frame of the computational model, and are prone to 
a long model construction time if the model is subjected to many driving inputs. 
The latter in particular is true in unsteady Computational Fluid Dynamic (CFD) 
applications where the moving solid boundary is often described by many structural 
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mode inputs. For example, a typical Boeing commercial airplane is modeled by as 
many as 200 structural modes. 

Recently, the Eigensystem Realization Algorithm (ERA)[1] was used successfully 
in application of CFL3D for aeroelastic flutter predictions [2], [3]. This method, which 
is usually used as a system identification technique, has a very attractive feature in 
that unlike model reduction methods based on Galerkin scheme [4] there is no need 
for on-line implementation of the algorithm. That is, it is a post-processing tool that 
identifies and generates system matrices based on the input and output data alone. 
Unfortunately, if the unsteady CFD model is driven by multiple structural inputs the 
computation time required to obtain all the pulse responses increases proportional to 
the number of the inputs, making the ERA slow and inefficient. 

Kim et al [5], [6] introduced the concept of the Single-Composite- Input (SCI) in 
applications of Karhunen-Loeve method to unsteady CFD models in time and fre- 
quency domains. The idea herein was that for a linear system, one can apply the 
multiple inputs simultaneously and get all the system responses that are necessary 
for the model reduction. Since the computational model needed to be executed only 
for the representative input, the model construction time was significantly reduced. 
In this paper, the same approach is adopted for fast and efficient model reduction 
of linear, finite-dimensional, discrete-time systems. To accommodate the SCI within 
the framework of ERA, it is necessary to modify the original algorithm. In particu- 
lar, it will be shown that the new formulation does not rely on the system Markov 
parameters explicitly. Instead, it performs the singular-value-decomposition (SVD) 
directly on the output measurements that are in general not attributed to pulse in- 
puts. Statistically independent random numbers are used in lieu of the pulses for 
the multiple input signals. Naturally, the new algorithm can also be used towards 
system identification provided that all the time measurements are available from ex- 
periments. Application of the SCI/ERA to computational fluid dynamic systems and 
formulation of reduced-order aeroelastic models are also presented. It is shown that 
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depending on how the displacement and velocity inputs on the moving boundary 
are treated, two different kinds of reduced-order aerodynamic and aeroelastic models 
could be generated. 

For a demonstration of the proposed method, two CFD models are considered. 
They are the Vortex Lattice Model (VLM) for inviscid, subsonic, incompressible flow, 
and the CFL3D for viscous, transonic, compressible flow. Reduced-order aeroelastic 
models are also constructed by combining the reduced aerodynamic models and the 
structural system. It is shown that not only the new method shortens the model con- 
struction time substantially, but the accuracy of the resulting reduced-order models is 
excellent. The proposed scheme has a great potential as a linear system identification 
and model reduction tool for large-scaled systems subjected to multiple inputs. 

2 Review of Pulse/ERA 

In this section, the Pulse/ERA (also known as ERA) is reviewed. For simplicity, 
only its fundamental state-space realization theory which is attributed to Ho and 
Kalman[7] is discussed. For a general description of the method, see Ref.[8]. It is 
assumed that the system under consideration can be described in the following finite- 
dimensional, discrete-time form: 

x" +1 =Ax" + Bu" (1) 

y" =Cx" + Du" (2) 

where 

x = (L x 1) state vector (3) 

u = (Ni x 1) input vector (4) 

y = (N 0 x 1) output vector (5) 

The system matrices (A, B), (A, C) are assumed controllable and observable. First, 
given M + 1 equally distributed time steps, t n = n At (n = 0, 1, 2, . . . , M), for a 
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single i — th input vector the system output is sampled subjected to the unit pulse, 
Assuming zero initial condition, x° = x(0) = 0, one obtains 

y? = d, 

y, 1 = cb, 

y- = CAbi 

y 3 = CA 2 b, 

yf = CA M1 bi (7) 

where 

hi = i-th column of B (8) 

dj = i-th column of D (9) 

The constant matrices in the above sequence are known as the Markov parameters[l]. 
This step is repeated Ni times for all inputs creating the sequence of pulse-response 
matrices: 

Y" ee [y? y 2 " ... y«J (n = 0,1,2,. ..M) (10) 

Next, based on the system Markov parameters define N Q x (Ni x (M - 1)) Hankel 
matrices: 

H 0 ee [Y 1 Y 2 ... Y^ 1 ] 

= C [I A A 2 ... A M " 2 ]B (11) 
Hi ee [Y 2 Y 3 . . . Y M ] 

= C [A A 2 ... A^B (12) 
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SVD of H 0 yields 

H 0 = UEV T 

= U fi E^E^VS (13) 
where R = rank(H 0 ). Finally, a balanced realization of the system under question is 
obtained by pseudo-inverting various submatrix components as 

D = Y° (N 0 x 1) (14) 

C ~ U fl £^ 2 (N 0 xR) (15) 

B ~ the first Ni columns of £^ /2 Y T R (R x JV,) (16) 

A ~ S^U^HiV^S^ 2 (RxR) (17) 

Since R « L, the above model represents a reduced-order realization of the original 
system. Note that the realization is optimal in that it is balanced between inputs and 
outputs. However, the total number of samples taken is AT* x (M + l), which increases 
proportional to the number of inputs. Also, for an accurate system identification H 0 
must have sufficient columns and rows, i.e., AT* x (M - 1) > R and N a > R. Figure 
1 shows a schematic of the Pulse/ERA procedure. 

For a very large data set Y" with many time steps and a large number of inputs, 
the Eigensystem Realization Algorithm/ Data Correlations (ERA/DC) method is pre- 
ferred to compress the amount of data and reduce the computation time required for 
the SVD of the Hankel matrix. See Ref.[l] for details. 

3 SCI/ERA 

The new method proceeds as follows. First, individual pulse responses are sampled 
for the first two time steps: 

Y° = [y? yS ... y° Ni ] (18) 
Y 1 = far} ... yy (19) 
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Next, construct the SCI as 



Ni 

^sci = E b* r ? (f° r states) (20) 

i=l 
Ni 

tfsci = E^ir" (for outputs) (21) 



r" = a sequence of arbitrary numbers (22) 

To ensure independency of the inputs, the signals must be as uncorrelated as possible. 
In an ideal case they would be statistically uncorrelated random signals, i.e., Cjj(m) = 
Elr™ r™~ m ] = 0 for i / j, but they are hard to construct for numerical analysis. 
Subject to the SCI, we sample the system response y n for n = 0, 1, 2, ... , M, and get 

y"o = Cx" 

= y"-Ey?rr (23) 

i=l 

= CAx" 

= y" +1 - Ey* 0 ^ 1 - Ey^r (24) 

Note that y™ 0 , y" x are measurements of the states in the reduced dimension of C and 
C A. Similar to the Hankel matrices, define 

H c0 = [yl 0 y c 2 0 ••• ycT 1 ] 

= C[x x x 2 ... x M - : ] (25) 

H c i = [y c \ y* x ... y^" 1 ] 

= CAfx 1 x 2 ... x M - : ] (26) 



SVD of H c0 yields 



H c0 ee USV T 

-p- ^[V s][3' 

= U B Si' 2 Ern (27) 
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Once again, R = rank(H c0 ). The size of the above matrices is N a x (M - 1), N t 
times smaller than the previous H 0 ,Hi. The new realization is then 

D = Y° (28) 

C ~ \J R J^ 2 (29) 

B ~ S^U^Y 1 (30) 

A ~ V R 1/2 V T R H cl V R V R 1/2 (31) 

Unlike the Pulse/ERA, the SCI/ERA is optimal in that it is balanced between states 
and outputs. As in the previous method, for an accurate realization H c0 must have 
(M - 1) > R and N a > R. However, the total no. of samples taken is only 
M + 1 + 2 x Ni which is much less than the previous iVj x (M + 1) when M samples 
of pulse response are taken for each input. Figure 2 shows a schematic of the SCI/ERA 
procedure. 

Two comments regarding the new algorithm are in order. First and foremost, 
compared to the Pulse/ERA the new algorithm requires a much smaller set of time 
measurements reducing significantly both the computation time and the bulk of the 
output data. Second, the new H c0 , H c i are constructed based on the sampled system 
states subjected to combined random inputs, and as such they are not directly related 
to the Markov parameters. However, the scheme does require the first two pulse 
responses yf and yj for each input in order to estimate the state meansurements 
according to Eqs. (23), (24). 

Although random signals are used exclusively in this paper, other types of signals 
can also be used for the SCI provided that they are statistically independent. It is 
well known that for a linear system any arbitrary response contains the fundamental 
characteristics under the assumption that the system is controllable and observable. 
It is this observation, along with the principle of superposition, that the present 
system identification scheme is based on. 



Boeing Proprietary 



10 



4 Alternative Scheme to Increase Measurements 

One of the requirements in the ERA methods is that for an accurate model con- 
struction one must have a sufficient number of output measurements, more than the 
number of singular modes that are extractable from the Hankel matrices. Failure to 
satisfy this requirement implies that we don't have enough modes to approximate 
the system response. When the requirement is not met, assuming again that (A, C) 
is observable we can expand the data matrices by sampling additional responses as 
follows. 



where 



H c0 i = 



H c n = 



C 
CA 
CA 2 

CA* 

y c i y c i 

y\K y 2 cK 
c 

CA 

CA 2 



[x 1 x 2 



yd- 1 



CA* 

Yc2 



yci 

V 2 

y C 2 



yd" 1 
yd" 1 



YcK+l YcK+1 • • • YcK 

y c \ = CA fc x" 

= y" +fc - E y°ir? +k - E yjrr*- 1 • 



(32) 



(33) 



Eyf^r 
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= y 

SVD of the new H c0 leads to 

H c01 ee U1E1V? 



EEyr 1 ^ 4 

j=i i=i 



(34) 



VI* 



from which we obtain 



the first N 0 rows of Y° 

the first N Q rows of Ui^S}^ 

s i a 2 u f R H c 1 1 Vi R 



where 



yi' 
y? +1 



y2 

y^ 1 



y^ 
y^T 1 



(n = 0,l) 



(35) 

(36) 
(37) 
(38) 
(39) 

(40) 



. yi y2 ••• Yn- 

The total number of measurements available is now (K + 1) x N 0 . It is noted that 
the additional time samples are required for the pulse response as well as for the 
response due to the SCI. More specifically, if K blocks of outputs are to be added 
pulse responses due to each input must be sampled at K additional time steps in 
addition the first two time steps, t = 0 and t = At. Also, the response due to the 
SCI must be sampled at K additional steps beyond the M-th step. The total number 
of samples to be taken is thus M + 1 + K + (2 + K) x Ni. K must be sufficiently large 
enough to satisfy the measurement requirement, (K + 1) x N 0 > R. Needless to say, 
this scheme requires extra computation time because of the additional time samples 
required in the data set. 
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5 Time Signals for Simultaneous Excitation 

A few candidate signals for the SCI are introduced in this section. Note that for an 
ideal linear model any of these SCIs can be used and all of them should yield ROMs 
of essentially the same quality. 

5.1 Random inputs based SCI (RSCI) 

It is natural to use random signals for construction of SCI. That is, use 

r" = a sequence of random numbers (41) 

in Eqs. (20) and (21). 

5.2 Filtered inputs based SCI (FSCI) 

To inject smooth inputs, one can filter the random signals through a low-pass filter. 
That is, 

r? = r]i 

= a sequence of filtered random numbers (42) 

Using the low frequency signals will allow better convergence when applying the SCI 
to CFD models. Furthermore, since the frequency content is limited it is possible to 
generate a smaller ROM directly from the SCI/ERA without any further reduction. 
A potential drawback is that if the filtered signals become too narrowly banded, they 
may not be as uncorrelated as desired. However, based on the theory of ergodicity [9], 
the statistical independence could be fortified by using longer signals and sampling 
the response for a longer period of time. 

5.3 Step inputs based SCI (SSCI) 

In analogy to a single step input, one can apply multiple step inputs in a sequential 
manner: 



Boeing Proprietary 



13 



= a step input applied at k r th step (43) 

where 

To assure independence of the inputs, the onsets of the signals must be apart from 
each other by a sufficient number of steps, i.e., k { must be large enough. 

5.4 Pulse inputs based SCI (PSCI) 

One can also apply multiple pulse inputs in a sequential manner: 

r? = &r ki 

= a step input applied at krth. step (45) 

where 

S n- ki = f 1 (n = ki) 

1 \ 0 (n = all other points) 

Again, A;, must be large enough to ensure independency of the applied signals. 

6 Second Reduction Based On FDKL/SCI Method 

In applications of discrete-time computational models, there exist two conflicting 
requirements for the incremental time step At. For numerical convergence one has 
to adopt a sufficiently small Ati. On the other hand, given the highest frequency 
of interest, u c , Nyquist criterion requires At 2 < Usually, for practical purposes 
At 2 » Aii. For instance, in a structural model that is coupled with a CFD 
model for aeroelastic applications, the highest mode usually has a much lower natural 
frequency than the highest frequency content in the aerodynamic model. If the signals 
used in the ERA methods are sharp as in the random, step, or pulse inputs, the 
SCI will excite all the system dynamics and hence this characteristic will be carried 



(46) 
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over to the ERA based reduced-order model. As a result, the ERA reduced-order 
model (ROM) is prone to have a large dimension to span the same high frequency 
range as the original full-order model (FOM), which suggests that there may be a 
room for further order reduction in the ROM. To perform a second order reduction, 
one can apply the Frequency-Domain Karhunen-Loeve (FDKL) method to the ERA 
ROM defined by matrices, (28)-(31), wherein frequency samples of the system within 
the given frequency range, (— uj c ,uj c ) are used to extract optimal modes, and a new 
reduced-order model is constructed via the Galerkin's approximation[10]. According 
to the method, the optimal or so called KL modes 0 4 are the eigenmodes of the 
covariance matrix K: 

= \i<J>i (47) 



K ti = Af(w«)*W r (48) 
ooi = sampling frequencies 

= [wiw 2 ... u M ] (49) 

where uj\ = — oo c and oo M = w c . X^ooi) are the frequency solutions of the ERA ROM 
subjected to the single-composite-input described by (20) and (21) except that it is 
prescribed in the frequency domain. Once the optimal modes are obtained, assume 

x ~ *p (50) 




(51) 
(52) 
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i?i is set to be equal to the rank of the covariance matrix which is usually smaller than 
R. After inserting (50) into (1) and (2) with the ERA ROM matrices, premultiplying 
both sides by * T yields a new reduced-order model as 

p n+1 = Ax p n + Bi u n (53) 
y" = Ci p n + D u n (54) 



Ai = * r A* (55) 
B x = $ T B (56) 
Ci = C* (57) 

The dimension of the new model is (Ri x Ri). 

7 Application to CFD Model 

In this section, application of the proposed model reduction/ system identification 
method to large-scaled CFD models is discussed. Unlike the general system described 
by Eq. (1) and (2), an unsteady fluid dynamic system is driven by displacement and 
velocity of its moving boundary surface simultaneously as they both contribute to the 
normal downwash on the surface. If one considers a statically nonlinear, dynamically 
linearized flow field, the unsteady fluid motion can be described as 

x" +1 = Ax" + B 0 u" + Biii" (58) 

y" = ?(Cx" + D 0 u" + Dili") (59) 

where 

x = (L x 1) fluid states (60) 

u = (JVj x 1) generalized displacements (61) 

li = (Ni x 1) generalized velocities (62) 
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y = (Ni x 1) generalized aerodynamic forces (63) 

q = dynamic pressure (64) 

It is noted that the above equations progress in non-dimensional time, r = rather 



than in the real time t and ( ' ) is the first derivative with respect to r. In fact, the 
dependency of the normal downwash on air speed disappears when the equations are 
discretized in r, as in Eqs. (58) and (59). The structural degrees of freedom, m are 
the generalized coordinates associated with structural modes. These modes are used 
to describe the motion of the lifting surface. Two different types of reduced-order 
fluid dynamic models can be obtained depending on how the inputs are treated. If 
necessary, the FDKL/SCI can be performed for additional reduction. 

7.1 Aerodynamic ROM with displacement and velocity in- 
puts 

One can treat u" and u" separately as independent inputs. Thus, for the pulse 
inputs 

< = < = '* ' {I (nJlitu)} W 
for i = 1,2,..., N i} we obtain the corresponding responses y°,y- at the first two 
time steps. Let us define 

y° = [y?y 2 °...y^y^ +1 y^ +2 ...y 2 Vj (66) 
y 1 = [y{yl..y 1 Nl y Ni+1 y Ni+ 2---y 1 2N i } (67) 

where the first Ni samples are due to the pulses in u" and the next N ones are due 
to the pulses in u". Next, we prepare the following inputs, 

Ni Ni 

h sci = E b oi^ + E h ur Ni+i (68) 

i=l i=l 
Ni Ni 

d 5c/ = £ do i»? + £ dur Ni+i (69) 

i=l i=l 
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Subject to the SCI, we sample the system response y" and get 

y"o = y" -Eyf»f ( 70 ) 

i=l 

y n c i = y n+1 - E y°< +1 - E yJr? (7i) 

i=l i=l 

Define 

h c0 = [yl 0 y c 2 0 ••• y^ 1 } (72) 
H cl ee [yh ... y^" 1 ] (73) 

where SVD of H c0 yields 

H c0 ee USV T 

<*4? !] [3] 

= U^E^EfvS (74) 
with ee ranA;(H c0 ). Hence, the reduced-order model is given by 

D 0 = the first JV< columns of Y° (75) 

Di = the second AT, columns of Y° (76) 

C ~ U fl Ef (77) 

B 0 ~ the first JV, columns of S H 1/2 U£ Y 1 (78) 

Bi ~ the second JV, columns ofS^ 1/2 U^Y 1 (79) 

A ~ E- 1/2 UiH cl V fl E- 1/2 (80) 

7.2 Aerodynamic ROM with displacement inputs 

The second scheme is to use only the displacements as the system inputs. This is 
possible by applying simultaneously the pulse and double pulse inputs, 
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and get the corresponding responses y^y^ at the first two time steps: 

yS = farSi y*2 ■■■ y° dNi ] (83) 

Y^ = [y^ yj 2 ... y^J (84) 

For a new SCI, we use 

Ni Ni 

*>sci = E b o^r + E b i^r (85) 
i=i i=i 

d n SCI = J^d^r? + Y,d u r? (86) 
j=i j=i 

where rf are discrete-time derivative of rf. To be consistent with the use of the 
double pulse defined in (82), r™ are obtained by filtering rf via <j", i.e., 

r? = conv(r!,8t) (87) 

which is equivalent to the backward difference equation, 



Subject to the new SCI we sample the system response y n and get 

y^co = y"-EySi»7 (89) 

i=l 

y d " c i = y" +1 - E yS^r +1 - E yi,r? (90) 
i=i i=i 

Defining 

Hdco = [yico ySeo ••• yfco 1 ] (9i) 
H dc i = [yLi ySei ••• yfcl 1 ] (92) 

the SVD of H dc0 yields 

H dc0 = USV T 

= U^sf V- (93) 
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with R = rank(Hdco)- The new reduced-order model has only AT* input channels 
and is in the form 

x" +1 = Ax" + Bu" (94) 

y" = ?(Cx n + Du") (95) 

where 

D = Y° (96) 

C ~ U*^ 2 (97) 

B ~ X£ 1/a U£Yj (98) 

A ~ V R 1/2 U T R H dcl V R V R V 2 (99) 

7.3 Reduced-order aeroelastic model 

We illustrate how aeroelastic systems can be constructed using the reduced-order 
aerodynamic models obtained in the previous section. 

First, we note that structural model is normally described in real, continuous-time: 

mii + cii + ku = y (100) 

( ' ) and ("), respectively, represent the first and the second derivatives with respect 

to t. Hence, to construct aeroelastic model the continuous-time euqation (100) is 
discretized in time: 

z" +1 = A s z n + B s y n (101) 

u" = C s z" (102) 

where 

Z - { I } ( 103 ) 

C s = [I 0] (104) 



Note that given At and V the consistent incremental time step, At = , must t 
used in the conversion to the discrete-time. 
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Aeroelastic Model I. 

In this approach, we treat u",ii n as independent inputs and apply the hg CI , dg CI 
given by Eqs. (68) and (69) at a reference dynamic pressure, q re f = \p re fV^ ef . The 
corresponding samples are taken and scaled by The ROM of the first kind de- 
scribed in Section 6.1 is then obtained by applying the SCI/ERA. An aeroelastic 
system that is valid at all V can be constructed by combining the resulting aerody- 
namic ROM with the structural equations: 

X n+1 = A dl X n (105) 

where 

X - { z } ( 106 ) 

Adl - [qB s C A s + qB s [D 0 D : ] J 
Denoting the eigenvalues of system matrix (107) as A d i», the aeroelastic eigenvalues 
in the Laplace domain are obtained as 

A CIj = ^ (108) 

For nutter instability, one must have Real(X cU ) > 0 and ||A d ij|| > 1 for any i. 

Aeroelastic Model II. 

One can also apply the ^sci^sci given by Eqs. (85), (86) and get the ROM of 
the second kind described in Section 6.2. This will produce aerodynamic system 
matrices, A,B,C, and D where B,D each has only AT* columns. The new resulting 
reduced-order aeroelastic model can be obtained as 

X" +1 = A d2 X" (109) 

where 

a - - [J.C a. + b 9 b:dc,1 ( iio > 
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Although q can change in Eq. (110), this model must be used only at the reference air 
speed V re f. It is well known that given a local angle of attack it, plunging rate ii, and 
local air speed V the total aerodynamic downwash at the moving boundary is ii + V u 
and as such it is impossible to separate out and account for the effect of V without 
having both the displacement and velocity channels. However, this drawback is easily 
remedied by adjusting the incremental time step according to At = when the 
air speed changes from one value to another and discretizing the structural model 
based on the new At. That is, if one leaves the V dependency in the structure, 



then the Aeroelastic Model II. becomes valid for all air speeds. 

8 Results and Discusion 

8.1 Unsteady vortex lattice example 

To demonstrate the proposed method, an unsteady vortex lattice aerodynamic 
model subjected to several structural mode inputs is considered. Both types of 
reduced-order aeroelastic models are constructed using the methods described in 
the previous section. For simplicity of implementation, only the results from the 
SCI/ERA method with a sufficient number of output measurements are presented 
here. Occasionally, the alternative scheme described in Section 4. was applied and 
indeed found to produce reduced-order models of essentially the same qualities as 
those presented here. 

Unsteady flow field around a flat, rectangular wing in incompressible, subsonic air 
flow is modeled by the vortex lattice formulation (Figure 3). The wing is 3 inches 
wide and 12 inches long, has 10 and 20 vortex elements in the chordwise and spanwise 
directions, respectively. The free wake has 40 and 20 vortex elements in the stream- 
wise and spanwise directions, creating a total of 800 degrees of freedom[13]. The 
wing structure is modeled using 6 vibrational (3 bending and 3 torsional) modes[14]. 




A BC S 

qB s (V)C A S (V) + qB s (V)BC s 



(111) 
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No structural damping was introduced at this time. Thus, the size of the full-order 
aeroelastic model is (812 x 812). 

For the Aeroelastic Model I. the reference air density and speed were set at 
1.23kg/m 3 ,80m/sec, respectively. The incremental time at this reference speed is 
At = ^ = 9.525 x 10" 5 sec. For the sampling of the vortex model, 480 extra 
outputs were extracted in addition to the 6 generalized aerodynamic forces at 481 
time steps. Applying 12 sets of random signal inputs simultaneously, 6 for u", 6 for 
u", yielded a single set of sampled data. 12 pulse inputs were also applied individ- 
ually at the first two time steps to generate Y° and Y 1 . Figure 4 shows three sets 
of the random numbers generated on Matlab. Out of 481 time samples, the SVD 
produced 413 linearly independent singular modes. This number was determined by 
the rank of H c0 matrix. Thus, the size of the aerodynamic ROM became (413 x 413). 
The reduced-order aerodynamic model was then coupled with the structural model 
to create (425 x 425) aeroelastic model (ROM I.). 

For the Aeroelastic Model II. the reference air density and speed were again set 
at 1.23kg/m 3 and SOm/sec. 6 sets of random signals and 6 sets of discrete-time 
derivatives of the random signals were applied for u" and li" using 481 time steps 
and the 486 output measurements. This yielded (329 x 329) aerodynamic ROM 
which when combined with the structural system, produced (341 x 341) aeroelastic 
model (ROM II.). It is noted that ROM II. is approximately 20 % smaller than ROM 
I. as a result of using only the half of the input channels. 

Next, the dimensions of the reduced-order aerodynamic models were further de- 
creased using the FDKL/SCI method. As mentioned earlier, the incremental time 
step embedded in both the FOM and SCI/ERA ROM is too small to be effective 
for various aeroelastic simulations which usually involve a low frequency range. Con- 
sidering that the highest free vibrational frequency of the structural modes is 4,160 
rad/ sec, the sampling range in the FDKL method was restricted to (—4, 500, 4, 500)rad/ sec. 
For the ROM I., out of 174 frequency samples within the range 129 KL modes were 
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selected based on the rank of the covariance matrix K. Hence, the size of the new 
reduced-order aerodynamic and aeroelastic models (ROM I.-FDKL) became 129 and 
141, respectively. Likewise, for the ROM II. 97 KL modes out of 130 frequency sam- 
ples in the same frequency range were selected yielding new (109 x 109) aeroelastic 
model (ROM II.-FDKL). For computational efficiency, these reduced-order models 
are to be preferred over the ROM I. and ROM II. 

Figure 5 presents (6 x 6) generalized aerodynamic forces obtained from the FOM, 
ROM I.-FDKL, and ROM. II-FDKL, in the nondimensional time, r. It is seen that 
despite the cut-off frequency range present in the latter two models, they reproduce 
the pulse aerodynamic responses of the original model very well. 

Figure 6 and 7 show two different scales of the aeroelastic eigenvalues of the various 
models in the s domain at V = 80m/ sec. It is seen that many eigenvalues of the 
reduced-order aeroelastic models match very well with those of the full model (Figure 
6). More specifically, the 12 complex eigenvalues associated with 6 structural modes 
agree very well between the FOM and ROMs, although the higher structural modes 
(5th and 6th) in the ROM II. and ROM II.-FDKL are slightly mismatched (Figure 7). 
Also noteworthy is that all the eigenvalues of the ROM I.-FDKL and ROM II.-FDKL 
are approximately within the specified bound, (-4, 500, 4, 500)rad/sec as the model 
was obtained based on frequency samples in the range. Presented in Figure 8 and 9 
are time responses of the first two structural modes due to an initial condition in tti. 
It can be seen that the three sets of curves are practically indistinguishable from each 
other. Next set of figures, Figure 10-13 show aeroelastic results at V = 121.2m/ sec. 
As can be seen from the figures, the wing is on the verge of flutter at this speed. Note 
how accurately the ROM I.-FDKL is able to reproduce neutrally stable, sinusoidal 
time responses (Figs. 12-13). However, the ROM II.-FDKL exhibits a noticeable but 
minor error in producing the transient response. 

Finally, the model construction time is compared between the two ERA methods. 
In order to obtain accurate and consistent singular modes, H c0 ,H c i matrices were 
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kept as square as possible by keeping the number of time samples approximately equal 
to the total number of measurements which is the sum of the number of generalized 
aerodynamic forces and the number of auxiliary measurements. The same numbers 
of time steps and auxiliary outputs were used both in the Pulse/ERA and SCI/ERA. 
Thus, in the first case where only the first bending mode alone excites the flow field, 
M = 131, iV 0 = 131. In the second case where the first bending and first torsional 
modes were included, M = 251, N 0 = 252, and in the third case where the first 
bending and torsional as well as the second bending modes excite the aerodynamic 
field, M = 281, N 0 = 283. 4 and 5 inputs were also used with M = 331, N 0 = 
334 and M = 411, N 0 = 415, respectively. Table 1. shows CPU seconds spent 
in constructing ROM I. on a SGI machine. Also presented in parenthesis are the 
dimensions of the corresponding reduced-order models. Table 2. shows CPU seconds 
consumed for ROM II. on the SGI machine. Presented again in Figure 14 is the 
CPU seconds vs. the number of inputs for ROM I. and ROM II. Note that these 
numbers represent total CPU seconds spent not only in sampling the response but 
also processing the data in the subsequent ERA schemes. As seen from the Tables 
and Figure, the new method clearly has an advantage over the Pulse/ERA in reducing 
the model construction time yielding saving factors of multiple numbers. Needless to 
say, as the number of inputs increases so does the saving. It is interesting that for 
a given number of inputs both ERA methods generate ROMs of very similar sizes. 
As expected, ROM II. of the SCI/ERA are as small as 80 % of the corresponding 
ROM I. Despite the different input channels, both SCI/ERA ROM I. and ROM II. 
require approximately the same CPU time implying that in the SCI/ERA the overall 
computation time is not very sensitive to the number of inputs. 

8.2 CFL3D example 

The SCI/ERA has also been applied to unsteady aerodynamic systems modeled by 
CFL3D code. CFL3D is a finite element program that based on the Navier-Stokes 
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equations models nonlinear viscous, compressible fluid motion in subsonic as well 
as transonic flow fields[15]. Although CFL3D describes statically and dynamically 
nonlinear flow, when subjected to a small amplitude it predicts dynamically linearized 
behavior of the flow field around the nonlinear static position in the form of Eqs. 
(58) and (59). For most practical aeroelastic analyses such as flutter prediction and 
dynamic gust loads calculation, the small amplitude approximation yields sufficiently 
accurate results. 

The aeroelastic formulation for the CFD code is slightly different than for the 
vortex lattice case in that the second aerodynamic model in section 6.2 is used but 
whenever the air speed V is changed the incremental time step At is adjusted ac- 
cordingly when discretizing the structural model, Eq. (100). As in the Aeroelastic 
Model I. the resulting model can account for the effect of changing the free stream 
speed. The airplane configuration under consideration is that of the Twin-Engine- 
Transport-Flutter-Model (TETFM). The aerodynamic grid is given by the so called 
the Wing-Pencil-Nacelle (WPN) model with the strut between the wing and nacelle 
omitted (Figure 15). The structural motion is described by 10 structural modes, re- 
sulting in total of 10 generalized aerodynamic forces per a mode shape. Figure 16 
shows the first 4 mode shapes. The computational aerodynamic model consists of 
approximately 700,000 cells and 30 blocks. For detailed description of the modeling, 
see Ref.[3]. Ref.[12] discusses the computational model and construction of the aero- 
dynamic and aeroelastic ROMs based on the SCI/ERA using various types of input 
signals mentioned earlier. In this section, only one set of ROM results that were 
obtained from the SSCI/ERA are presented from the cited reference. 

The WPN model was examined at M = .831. With At = 3.34 x 10" 4 sec and 
995 time steps, aerodynamic ROMs were created using both the Pulse/ERA and 
SSCI/ERA methods based on the displacement method described in section 6.2. The 
size of the aerodynamic ROMs is (512 x 512). Given the 10 mode inputs, the total 
number of time steps consumed in the time marching required for the SCI/ERA 
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process was 995 + 81 x 10 = 1805. On the other hand, the traditional ERA required 
total of 10 x 995 = 9,950. Using a single CPU of IBM/Regatta machine, the total 
CPU hours for the Pulse/ERA calculation was 366 while the SSCI/ERA required only 
64 hours. Thus, the computational cost was reduced by a factor of 5.7. Figure 17 
shows the V-g plot of the two different aeroelastic models. As seen from the figure, 
the two ROMs match very well in flutter characteristics predicting the two flutter 
points within just 0.13 % and 0.36 % differences, respectively. 

9 Concluding Remarks 

In this paper, an efficient time-domain model reduction/system identification tech- 
nique has been presented and demonstrated for linear dynamic systems that are 
subjected to multiple right-hand-side inputs. The method is 100 % a post-processing 
that does not require modifying the original code and takes only input and output 
data for the model construction. The new method is based on a direct singular value 
decomposition of the output measurements that are not necessarily attributed to 
pulse inputs but due to multiple signal inputs applied simultaneously at the input 
channels. Compared to the Pulse/ERA, the SCI/ERA can significantly reduce the 
model construction time and compress the amount of output data. Therefore, it is 
very attractive for large-scaled dynamic systems with multiple driving inputs such as 
CFD models wherein the moving boundary input is often described by many struc- 
tural modes. For practical applications a second model reduction of the ERA ROM is 
desirable. For this purpose, the FDKL/SCI method is recommended for its efficiency 
and accuracy. 

Acknowledgment 

The author is grateful to his manager, Steven R. Precup and Senior Technical Fellow, 
Kumar G. Bhatia for their continuous support for this research. Note that the algo- 
rithm presented in this paper is a Boeing Intelelctual Property under consideration 



Boeing Proprietary- 



Table 1: CPUs of ERAs Applied to VLM (ROM I.) 



Ni 


Pulse/ERA 


SCI/ERA 


1 


12.3 sec (92) 


6.7 sec (89) 


2 


74.8 sec (180) 


16.6 sec (182) 


3 


150.8 sec (221) 


20.6 sec (226) 


4 


327.8 sec (297) 


28.2 sec (298) 


5 


762.3 sec (336) 


43.7 sec (340) 


6 


1,525.5 sec (395) 


63.2 sec (413) 



NOTE: Number in ( ) is the size of ROM. 



Table 2: CPUs of ERAs Applied to VLM (ROM II.) 



N 


Pulse/ERA 


SCI/ERA 


1 


6.3 sec (92) 


6.8 sec (89) 


2 


32.8 sec (169) 


17.00 sec (167) 


3 


65.9 sec (215) 


21.0 sec (214) 


4 


130.1 sec (258) 


27.6 sec (257) 


5 


279.6 sec (304) 


44.1 sec (305) 


6 


540.4 sec (316) 


61.4 sec (329) 



NOTE: Number in ( ) is the size of ROM. 
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Figure 1: Schematic of Pulse/ERA process. 
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Figure 2: Schematic of SCI/ERA process. 
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Figure 3: Vortex lattice grids for rectangular semi-span. 
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Figure 4: Statistically independent random signals. 
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Figure 5: Generalized aerodynamic forces at V 
I.-FDKL (129) ... ROM II.-FDKL (97)). 
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Figure 6: Aeroelastic eigenvalues at V = 80m/sec. 
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Figure 7: Aeroelastic eigenvalues at V = 80m/ sec. 
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Figure 8: Mode 1 aeroelastic response due to IC in mode 1 at V = 80m/ sec. 
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Figure 10: Aeroelastic eigenvalues at V = 121.2m/sec. 



Boeing Proprietary 



40 



4000 
3000 
2000 
1000- 

| <>-♦! 

-1000- 
-2000 - 
-3000 
-4000 



o FOM(812) 

ROM I. (425) 
□ ROM II. (341) 
<> ROM l.-FDKL (141) 
a ROM ll.-FDKL (109) 



-100 -80 

Figure 11: Aeroelastic eigenvalues at V = 121.2m/ sec. 
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Figure 14: Model construction time of VLM ERA ROMs vs no. of inputs. 
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Figure 15: The CFL3D grids for WPN model. 
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Figure 16: First 4 structural mode shapes for WPN model. 
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Figure 17: V-g plot of WPN model: Pulse/ERA vs. SSCI/ERA (growth rate only). 



